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Abstract 

Physicists have used billiards to understand and explore both classical and quantum chaos. Re- 
cently, in 2001, a group at the University of Texas introduced an experimental set up for modeling 
the wedge billiard geometry called optical billiard in two dimensions. It is worth mentioning that 
this experiment is more closely related with classical rather than quantum chaos. The motivation 
for the present work was born from the idea of laying the foundations of a quantum treatment 
for optical billiards, named "The Escape Problem", by presenting the concept of a Transparent 
Boundary Condition. We consider a "gas of particles" initially confined to a one dimensional box of 
length L, that are permitted to escape. We find the solution of a Quantum Initial Value Problem 
using a numerical method developed and entirely checked with an exact, analytic method. The 
numerical method introduces a novel way to solve a Diffusion Type Equation by implementing 
Discrete Transparent Boundary Conditions recently developed by mathematicians. 

PACS numbers: 
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I. INTRODUCTION 



The seminal model for understanding both classical and quantum chaos is the billiard. 
Depending on different considerations, such as shape and topology, a billiard model can 
present either stable or chaotic behavior. The most celebrated versions of classical billiards 
in two dimensions, in which chaos has been observed, are the Sinai billiard^, Bunimovich 
stadium^, and Polygonal billiards'^. In all of these examples, the shape of the boundaries 
plays a crucial role. In recent years, physicists have detected chaotic behavior at the nano 
scaled making an important impact in technology and, once again, billiards have helped 
to model quantum devices^ such as nanotubes, quantum dots, etc. which are governed by 
quantum theory. 

A novel model of billiard is the one introduced in 1986 by Lehtihet and Miller^ known 
as the "Wedge Billiard." The model consists of a symmetrically inclined wedge of angle 
29 with respect to the direction of a constant gravitational field g in which the particle 
is confined. They found surprising properties when the parameter 6 is changed. Later, 
in 2001, Valery Milner-, working in Mark Raizen's laboratory at the University of Texas, 
introduced an experiment referencing the Wedge Billiard model geometry and confirming the 
properties mentioned above. This first experiment of billiards was named " Optical Billiards.'" 
It is worth mentioning that this experiment is more closely related to classical rather than 
"quantum" chaos. While very low by normal standards, the temperature they employed is 
not low enough to easily exhibit any quantum effects. Consequently, a classical theoretical 
model was adequate to obtain a good point of comparison. 

A key element in the experiment is that, by removing a small segment of the boundary, 
an escape route is provided at the billiard vertex. As time progresses, the atoms, moving 
under the influence of gravity, will eventually exit the billiard. A measure of the influence 
of "chaotic" orbits is provided by the mean lifetime that the atoms remain in the billiard. 
This scenario is mimicked in this work by setting up the Escape Problem {EP). Currently 
experimentalists are starting to probe the quantum regime. 

The idea of the work presented here springs from laying the foundation of a quantum 
treatment for optical billiards. However, the quantum mechanical problem is much more 
difficult since the system is two-dimensional (equivalent to four dimensions in the phase 
space), and not integrable. Thus analytical methods are not available, and a viable numerical 
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method for solving the Schrodinger equation is required. Since it is extremely difficult 
to develop a numerical method for a two dimensional integro-differential equation with 
nonlocal transparent boundary conditions and skew boundaries, this work will develop an 
essential numerical method for obtaining the solution of the non-relativistic one-dimensional 
Schrodinger equation. 

The most demanding part of this work is the implementation of Transparent Boundary 
Conditions (TBCs)^. The TBCs, recently developed by mathematicians, arise in the ne- 
cessity to deal numerically with the natural infinite domain of the wave function. That is, 
due to the limitation of the computer core size, the infinite wave solution of the Schrodinger 
equation has to be solved in a finite domain by imposing artificial boundary conditions. The 
criterion of transparency is that the incident wave at the boundary has the smallest reflection 
coefficient possible. If the solution with these artificial conditions agrees with the infinite 
solution, the artificial boundary conditions are said to be transparent. From the several 
available approaches to derive TBCs, this work is based on the excellent treatment carried 
out by Anton Arnold^and his student Mathias Earhardt^*^ that concerns the transport of 
a quantum particle that enters one side of a finite domain and exits from the opposite side. 

In the following, in the next section, we ffist set up the "Escape Problem, i.e. the escape 
of a particle from a finite region, as a Quantum Initial Value Problem (QIVP). We then 
develop two different approaches for solving the EP. The ffist is the analytic method that 
provides the certainty of the result. The second is the numerical method, where an algorithm 
will be meticulously presented by introducing Discrete Transparent Boundary Conditions 
(DTBC)^, based on the implicit Crank- Nichols on method. Subsequently, the consistency 
of both methods is conffimed in section III by comparing the solutions for the real and 
imaginary parts of the wave function at different times, showing an excellent agreement. In 
the last section, conclusions are drawn from the work that has been accomplished, stressing 
the valuable physical information provided by the wave function as well as the numerical 
importance of the method developed in this work. 
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II. THEORY 



A. DESCRIPTION OF THE SYSTEM. 



The system representing the one-dimensional Escape Problem [EP] is a group of particles 
restricted inside a region of size L delimited by boundaries and free of external influences 
where, abruptly, one of the boundaries becomes transparent, allowing particles to escape. 
The group of particles is a dilute gas rarefied enough to be considered as a "Knudsen Gas"— 
whose density is so small that only the interactions with the boundaries are relevant. From 
the perspective of quantum mechanics, the EP is a QIVP. Thus, the intention is to solve the 
1-D Schrodinger equation from a given initial condition. The EP is similar to a diffusion type 
problem; it has the same structure as propagation problems in which a partial differential 
equation is solved with the aid of a known initial value. The evolution of the gas inside the 
region < x < L will be determined by the wavefunction at the initial time. To demonstrate 
the approach, here the initial wavefunction is chosen to be the ground state of an Infinite 
Well Potential of length L. Therefore, the initial value or the initial condition inside of the 
box is 

■qfj{x,0) = sm{kx) (1) 

with k = . Of course, outside of the box the value of the initial wave function vanishes. 
The evolution of the EP is represented in the figure [H 



B. ANALYTICAL METHOD. 



The non-relativistic, time dependent, one-dimensional Schrodinger equation is 

dt 2m dx"^ 

with the initial condition: 

^(x,t)|t=o = ^/(a;,0). (3) 

The Laplace Transform Method is one of the many appropriate techniques to solve differ- 
ential equations. The first task is to determine the Laplace transform of the wave function 
and its time derivative. For this particular example the definition of the Laplace transform 
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is: 



oo 



x,s) = e 







or, 



{x,s) = C{<ilix,t)}. (4) 



C\^^^^} = sct>{x,s)-^j{x,0). (5) 



By integrating by parts, 



dt 

After using equations (jlj) and ([5]), the Schrodinger equation ([2]), after some algebraic 
re-arrangements, becomes: 

^-^^^^^-^ + ias0 (x, s) = zq;\1// (x, 0) (6) 

with 

A Green's Function can be used to solve the inhomogeneous differential equations (E]). 
The following notation and definitions are necessary before using the Green's Function. 
The Laplace transform of the Green's function is denoted as: 

C{G{x,x',t)} = g{x,x',s) (8) 

where G{x,x',t) is the solution in the non-transformed space. The Green's function for 
equation ([6]) satisfies 

— ^ ^ ' ' ^ + iasg (x, x' , s) = 6 (x — x') (9) 
dx'' 

and its solution is given by 

POO 

(f){x,s)=ia / dx'g {x,x' , s)'^j {x' ,0) . (10) 
Jo 
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The complete time dependent solution of the initial value problem \1/ {x, t) is obtained by 
inverting the Laplace transform: 

/•CO 

m{x,t) = iaj dx'G{x,x',t)^i{x',0) (11) 
Jo 

The remaining procedure consists in constructing the Green's function g{x,x',s) for our 
particular problem, establishing G{x,x',t) , and finally substituting it into equation (fTTj) . 

In order to construct the Green's Function g{x,x',s) , beginning from equation ([9]) for 
x ^ x' and assuming ReP > , the form of the postulated solution is g = e^^^. Thus, 




For X < x' 



And for x > x' 



P=J^il-t). (12) 



x = (0) = 

3 + = (13) 

91 = i?(e^--e-^-) (14) 

92 = Ae-""^ (15) 



The following two conditions give A and B. For continuity gi = g2, and for the disconti- 
nuity of the derivative ^L.' - ^L' = 1- 

•' ax '-^ ax '-^ 

These two conditions lead to the value of: 

B = - — e-^"'. (16) 
Using this value of B, the value for A is derived as 

A = ^ (e-' - e--') . (17) 
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The equation (fT^ and the equation (fT5|) can be rewritten to obtain: 



Alternatively, 



and. 



for X < x' 

9i = 
for X > x' 



2P 

1 



= -^(e^^'-e-^^')e-^^ (19) 



9ia = -^e-^(^'-^) (20) 
91b = T^e-^^-^-') (21) 



= -^e-^(-^') (22) 
^2. = 7^e-^(^+^'). (23) 



In general, the inversion of the Laplace transform of any of the g functions has to be done 
via contour integration in the complex s plane. In addition, extensive mathematical tables 
are also available with many worked out examples^. This section shows only one of the g 
functions, gia , as an illustration of the use of these tables (the rest of functions gib, g2a, g2b 
are easy to obtain due to their similarity). The following expressions can be found^ 



7i' 



Fit) 
1 
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for A; > 0. Then, using ([?[) and (fT2|) it is easy to see that 



9la 



2fl-i) 



I ms 

h 



(24) 



and its corresponding inverse Laplace transform 



G 



= -3t[x(l-*)'(^'+^)] 



lb 



2(l-z) Artf 



(25) 



Similarly, it is easy to obtain 



G 



2a 



(26) 



G. 
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e-s[f(i-^)'(-'+-)] 
2(l-z)^Atf 



(27) 



Defining the following constants C 



2(l-i) 



and a = ^ the total Green's function is 



found to be: 



G (x.a;', t) 



(28) 



Once the Green's function is constructed by the above procedure, in order to obtain the 
wave function it is necessary to substitute equation (l28l) into the equation (fTTl) resulting in 
the complete solution: 

m t) = z^ dx' (-eT(-'-)^ + vl>, (x', 0) 



Mathematica was employed for these computations and for plotting the results. 



C. NUMERICAL METHOD. 

As mentioned earlier, the QIVP can be treated as a propagation problem. These are 
initial-value problems governed by a parabolic Partial Differential Equation {PDE) of first 
order in the time. A familiar representation of a parabolic PDE in one dimension is the 
Diffusion Equation, 

ft = af,,. (29) 
8 



Subscripts t and xx represent the first time derivative and second position derivative of the 
function /. The factor a in the equation is called the diffusivity which is defined by the 
system under investigation. The Schrodinger equation (in dimensionless units) using the 
same notation as in equation f l2^ becomes 

m>^ = -]^m^^. (30) 



Except for the imaginary number i, the equation (130|) is identical to fl29|) . Therefore, it 
is mathematically correct to proceed to solve ( !30l) with the complex extension of the same 
tools used in the numerical method for the Diffusion Equation. Numerical methods^^ solve 
the PDF by transforming the integral problem into an algebraic one that is computationally 
accessible. The Crank- Nicholson method is the preferred numerical algorithm used to solve 
the Schrodinger equation as a diffusion type equation. It is an implicit algorithm valid 
through second order in both space and time coordinates, so it is very stable^. 

To make that equation equivalent to the escape problem, it is necessary to employ the 
mathematical concepts of TBCs. The discretized domain of the solution is constructed using 
the following finite difference grids 

Xi = {i- 1) Ax (31) 
r = n/\t {At constant) (32) 

with i = 0, 1, 2, . . . , imax and n = 1,2, ... , nmax . The domain D {x, t) is from to L on 
the x — axis and the solution is marching in a positive time direction on the y — axis. We 
introduce the following notation in the context of the finite difference method grid point 
{i, t) — ^ {xi, t^) , function / (x^, t^) — )■ f^, first time derivative — )■ ft\2, and second space 

derivative > fxx\l:- 

The 2nd Order Central Space approximation of the second derivative is: 

in o £n I £n 

Jxxli - {'i'i) 

To obtain second order precision in the time, a key point of the method is to estimate 
derivatives at half integral time steps. The 2nd Order Central Time approximation of the 
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first derivative is given by: 

ftC = (34) 

The second space derivative is 

fxx\i ~ = 2 ^fxx\i "I" fxx\i ) • (35) 

The future time level occurs at n, the past time level at n—1 and the algebraic relationship 
for the one-dimensional approximation of the Schrodinger equation is: 



Then the Crank- Nicholson difference equation is: 

- ^7+1 + C*^ - = ^'I'l - C^r^ + ^tl (36) 



with C = (2 - ip) , C = (2 + zp) , and p 



The algebraic relationship develops a tri-diagonal system represented by the following 
equation 

A^" = 6, (37) 

of a set of imax simultaneous linear equations, where imax is the maximum number of grid 
points , imax — 1 equations come from the interior scheme and one equation comes from the 
Right Transparent Boundary Condition (at i = imax ). The left boundary at i = is still a 
Dirichlet Boundary Condition and it is set to zero in order to force the function to stay in 
the domain. 

The code presented in this work was based on the design of DTBCs derived by the work of 
Arnold^i^ with the inclusion of modifications that adapt DTBCs to the particular example 
of the escape problem. The original problem as conceived by Arnold and derived in the 
dissertation of his student, Earhardt^^, evaluates the transport of a quantum particle that 
enters one side of a finite domain and exits from the opposite side. 

The Right TBC is: 
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In contrast with Dirichlet conditions, they are not stagnant in time and the numerical 
method requires the storage of the past history values at all time levels at the boundaries. 
The convolution term on the right hand side of the equation (|38|) appears as a fractional (|) 
time derivative. Derivation of the DTBCs uses the Crank-Nicholson scheme at a discrete 
level. The discrete TBCs for the one dimensional Schrodinger equation for n ^ 1 on the 
right is at j = J 

n-l 

fc=i 

with, 

1- = (1 _ ,| + ^) ^0 + (1 + ,| + ^) ^1 + ^^^-in^,MH)_f^^^ 



^ - h~Kr^ ^1 - arctan 7^234^—^, 

3 

LLj = —r, r=, = f , CTj = 2 Ax Vj , 



J 2Y V-''' ' 

The Right discrete boundary condition is derived at i = imax — 1 

_ Vl/" + C\lf" 1 - n = (39) 

*rJ. - C'vl>r-i_i + ^La._2 = (^max - 1) (40) 

where 

z = L, vI/«_. = v[; (L, t) , (41) 

the Right Transparent Boundary Condition, is part of the system of equations. The Discrete 
Transparent Boundary Condition (138|) is: 



imax—1 imax imax 



^,n+l _ ,(0) ,T,„+1 

n-l 



E ^taMmax ' ^La.-i = {^rnax) . (42) 
fc=l 

Note that {imax) will change at every time level n. This is the value at that boundary 
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that must be stored by the selected numerical method. For the purpose of illustrating the 
procedure the first n — 3 equations are: 



imax—l 
imax—1 



imax—1 



imax imax imax—1 



;2 ^1 +1^ 

imax imax ' "rmirr- ^ 



imax imax 



^3 

imax—1 

imax— 1 



n = 1 

imax imax 

— 6° (imax) 
n = 2 

(imax) 
n — 3 



imax imax 



— (imax) 



(43) 



(44) 



(45) 



In matrix form, 



" c 


-1 

















-1 


c 


-1 

















-1 


c 


-1 



















-1 


c 


-1 




imaa;— 1 


_ 










-1 


/(O) 
imax 







6(2) 

6(3) 
6(4) 

6 {imax — 1) 
6 (imax) 



(46) 



A code called RTBC was written in C which employed the LU Factorization package to 
solve the complex valued tri-diagonal system of simultaneous linear equations numerically. 
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III. RESULTS. 

Since we use the implicit method called Crank- Nichols on with a degree of approximation 
through second order that is unconditionally stable, there is no formal restriction on the 
step sizes, and we expected little difference between the results obtained by both methods. 
The real and imaginary parts of the wave function are plotted for a set of different values of 
the time. For simplicity, dimensionless units are used where m = h = 1 . In these units we 
choose the length of the region to be L = 2. The set begins with two small values of time 
and ends with a relatively longer value compared with the value of the natural period. In 
dimensionless units, T = — = 5.092, for the first ground state with a value of a; = 1.233. 

A. ANALYTICAL METHOD. 

The exact, time dependent wave function on the semi- infinite domain is obtained by 
replacing the value of the initial condition in the equation f l29p as follows: 



After performing the integration using Mathematica Software, the explicit form of the 
wave function is 




(47) 




{Erfi 




—e L Erfi 



(i + i) {irt + aL {L - x)) 
y/aL\/i 



—e ^ hrji 



(i + I) {-TTt + aL {L + x)) 




+ Erfi 




)] 



(48) 



where Erfi is the error function with imaginary argumenlr^. 
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B. NUMERICAL METHOD. 



A code based on the Crank- Nichols on algorithm was developed to solve the TBC problem. 
All of the code was written in the C programming language and the calculations were 
performed on a LINUX workstation using the GNU compiler. The major effort has gone 
into adapting the techniques of TBCs because it is complicated by the fact that the boundary 
condition is non-local in the time, but rather depends on a convolution related to fractional 
calculus that cannot be treated by standard techniques. The standard algorithm had to be 
modified to include the non-local, time dependent boundary conditions, which now have to 
be updated at each time step. To solve the Schrodinger equation on the finite interval, the 
LU decomposition method was employed. The function libraries were obtained from the 
open source literature for the GNU compiler. 

C. COMPARISON OF THE METHODS. 

The first set of plots compare the real part of the wave function for both methods at the 
same time steps t = 0.1 , t = 0.2 , and t = 1.4 as shown in figures |2] and [31 

The second set of plots, figures H] and |5l compare the imaginary part of the wave function 
for both methods at the same time steps t = 0.01 , t = 0.1 , and t = 1.4. 

The comparison between the plots for each method shows an excellent agreement as we 
expected. This confirms the adequate implementation of the numerical method for a QIVP. 
Therefore, the basis can be established for a more complex version of the RTBC. 

IV. CONCLUSIONS. 

In this work we investigated the problem of solving the time dependent Schrodinger 
equation within a confined region when a constraint is suddenly removed. We accomplished 
two important goals; the first one is of physical significance. We determined the wave 
function at any time t, and therefore provided all the physical information for the evolution 
of the EP within the trap. Therefore, we are able to extract important physical properties 
of the system such as the mean fraction of mass inside the region, which is called the survival 
probability, as well as the time dependence of the properties of the particle kinetic energy, 
potential energy, the total energy, etc.. Since it is assumed that the particle density in the 
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system is small, interactions can be ignored (it is a Knudsen Gas-^, and multiplying the 
survival probability by the total number of particles gives us a quantitative value of mass 
inside the region. 

The second conclusion concerns the numerical significance. We have shown that the 
discretization developed by Arnold and Ehrhardt for the solution of open boundary prob- 
lems provides a viable numerical approach for solving the time dependent Schrodinger 
equatioii^""— when a constraint is suddenly removed. Although the present work consid- 
ers a one-dimensional model that is not sufficient to exhibit chaos, the code design forms the 
essential basis of future research on non-integrable, two-dimensional billiard models where 
chaos is present. The influence of chaos in these models will be explored by measuring the 
survival probability of particles in a trap with an open boundary. 

The results obtained with the code presented here show an excellent agreement with the 
analytic approach, offering reliability not only in the result obtained but also in the novel 
numerical method used. The success in the development of the numerical method opens 
up a reasonable extension to higher dimensional models. Also, it offers the possibility of 
understanding the quantum mechanical version of billiard models that experience classical 
chaotic behavior. In spite of the very low, by normal standards, temperature regime, a 
classical theoretical model was still adequate for the analysis of the Austin experiment. The 
quantum regime should also be accessible experimentally at yet lower, but still attainable, 
temperatures. Experimentalists are beginning to probe this regime. 
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FIG. 3: Real part of the wave function obtained numerically. 
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FIG. 4: Imaginary part of the wave function obtained analytically. 
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FIG. 5: Imaginary part of the wave function obtained numerically. 
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